#ifndef BOUNDCONDITIONS_H
#define BOUNDCONDITIONS_H

#include <functional>
#include "Vec.h"
#include "FunctionSet.h"

template <int N>
using Point = Vec<Real, N>;

using Point2d = Vec<Real, 2>;

template <class T>
inline int sign(const T &x)
{
    return (x > 0) ? 1 : ((x < 0) ? -1 : 0);
}

template <int N>
struct SingleBoundaryCondition
{
    /**
     * @brief \phi(x) = 0 represents the boundary itself.
     *      \phi(x) > 0 represents the interior of a region.
     *      \phi(x) < 0 represents the exterior of a region.
     * 
     */
    std::function<Real(Point<N>)> phi;

    // \alpha * u + \beta * \par u = g;
    double a;
    double b;
    std::function<Real(Point<N>)> g;

public:
    SingleBoundaryCondition(std::function<Real(Point<N>)> phi, std::function<Real(Point<N>)> g, double a, double b) : phi(phi), g(g), a(a), b(b) {}
    
    /**
     * @brief Tell the location of a point w.r.t the boundary.
     * 
     * @param p 
     * @return int : 0 on the boundary , 1 inside the region , -1 outside the region.
     */
    int Location(const Point<N> &p) const
    {
        return sign(phi(p));
    }

    /**
     * @brief Compute the partial differentiation 
     * along one coordinate axis of the boundary function \phi.
     * 
     * @param p 
     * @param axis : the axis along which partial diff is computed
     * @return Real : the value of partial diff
     */
    Real partial(const Point<N> &p, int axis) const;

    /**
     * @brief the outer normal vector on the boundary point p
     * 
     * @param p 
     * @return Vec<Real, N>
     */
    Vec<Real, N> Normal_Vec(const Point<N> &p) const;
};

template <int N>
struct BoundaryConditions
{
    //all single boundarys form the whole closed boundary
    std::vector<SingleBoundaryCondition<N>> bounds;

public:

    BoundaryConditions() = default;

    BoundaryConditions(std::vector<SingleBoundaryCondition<N>> &&bounds) : bounds(bounds) {}

    BoundaryConditions(const FunctionSet<N>& funcs,double a,double b)
    {
        for (int i = 0; i < 2 * N;++i)
            bounds.push_back(SingleBoundaryCondition<N>(
                funcs.Boundary[i], Function<N>{},
                a, b));
    }


    /**
     * @brief Tell the location of a point w.r.t the boundary.
     * 
     * @param p 
     * @return int : 0 on the boundary , 1 inside the region , -1 outside the region.
     */
    int Location(const Point<N> &p) const
    {
        for (auto &bc : bounds)
        {
            int location = bc.Location(p);
            if (location <= 0)
                return location;
        }
        return 1;
    }
};

template <int N>
Real SingleBoundaryCondition<N>::partial(const Point<N> &p, int axis) const
{
    double dx = 1e-8;
    Point<N> p1 = p, p2 = p;
    p1[axis] += dx;
    p2[axis] -= dx;
    return (phi(p1) - phi(p2)) / (2 * dx);
}

template <int N>
Vec<Real, N> SingleBoundaryCondition<N>::Normal_Vec(const Point<N> &p) const
{
    Vec<Real, N> nor_vec;
    for (int i = 0; i < N; ++i)
    {
        nor_vec[i] = -partial(p, i);
    }
    return normalize(nor_vec);
}

#endif